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We study a system of trapped bosonic particles interacting 
by model harmonic forces. Our model allows for detailed ex- 
amination of the notion of an order parameter (a condensate 
wave function) . By decomposing a single particle density ma- 
trix into coherent eigenmodes we study an effect of interaction 
on the condensate. We show that sufficiently strong interac- 
tions cause that the condensate disappears even if the whole 
system is in its lowest energy state. In the second part of our 
paper we discuss the validity of the Bogoliubov approxima- 
tion by comparing its predictions with results inferred from 
the exactly soluble model. In particular we examine an energy 
spectrum, occupation, and fluctuations of the condensate. We 
conclude that Bogoliubov approach gives quite accurate de- 
scription of the system in the limit of weak interactions. 

PACS number(s): 03.75.Fi, 05.30. Fk 



I. INTRODUCTION 

Recent advances in the trapping techniques have re- 
newed interest in various aspects of many body theory. 
In fact a cloud of weakly interacting trapped atoms is 
an ideal system for which various aspects of many body 
theory can be tested and verified. The ideal bosonic gas 
undergoes the Bose-Einstein condensation if the phase- 
space density exceeds one. This phenomenon manifests 
itself by the macroscopic occupation of the single par- 
ticle ground state. In the case of an interacting system 
the condensate wave function can be defined by the spec- 
tral decomposition of the one-body density matrix. This 
decomposition is closely related to the off-diagonal long 
range order Jl| or the existence of the order parameter, 
i.e., the 'classical' field with given amplitude and phase 
commonly used in the theory of superfluidity |^|J3| . Real- 
ization of decomposition procedure is practically impos- 
sible because it requires a full solution of the many-body 
problem. Mean-field approaches are commonly used in- 
stead. The basic idea for a mean-field description of the 
dilute, weakly interacting Bose gas below transition tem- 
perature was introduced by Bogoliubov jf|. Most of the 
results for the interacting Bose-Einstein condensate are 
obtained within the Bogoliubov theory which in many 
cases provides a reliable quantitative description of the 
quantum Bose gas. Indeed, the low energy excitation 
spectrum of the trapped condensate well below transi- 
tion temperature Q agrees remarkably well with predic- 
tions based on the Bogoliubov theory ||. On the other 
hand the mean-field Bogoliubov approach fails to repro- 



duce excitation spectrum at higher temperatures - close 
to the transition point 0. Therefore, the question about 
limits of validity of the Bogoliubov method is of great 
importance. One possible way to test the quality of this 
approximation is to go beyond the mean field theory. An- 
other possibility to assess the usefulness of Bogoliubov's 
theory is to study the exactly soluble models and to com- 
pare their predictions with those based on the approxi- 
mate method. The models provide not only a unique 
soluble many body problem but also allow to verify var- 
ious approximations. This will shed some light on the 
validity and exactness of the Bogoliubov method, widely 
used in many body physics ||. 

There are only few exactly soluble models of quantum 
systems where the interactions between atoms is cho- 
sen in the form allowing for the exact analytic solution. 
These are: (i) the one-dimensional model of impenetra- 
ble bosons introduced by Girardeau 0, (ii) its contact 
potential version formulated by Lieb pQ | ; (iii) the model 
of particles interacting by harmonic forces p~j[— |T^] . Al- 
though in the first two cases the formal solution is given 
but in practice the problem is still quite complicated and 
quantitative calculations can be done for a very small 
number of particles only [Q. The latter case seems to 
be much simpler because, as it has been shown in [ jl5| , 
it can be reduced to the problem of noninteracting par- 
ticles in a harmonic trap. Therefore in the following we 
are going to examine, within this exactly soluble model, 
various concepts and methods related to the interacting 
Bose-Einstein condensate. 

This paper is organized as follows. In Sec. II we present 
exact results regarding properties of harmonically inter- 
acting bosons trapped within harmonic potential. The 
results are obtained within a soluble model that was de- 
veloped in p5| . In the second part of Sec. II we find the 
analytic expression for the order parameter and study the 
effect of quantum depletion of the condensate as well as 
quantum fluctuations at zero temperature. In the third 
part we analyze the thermal properties of the system. In 
Sec. Ill we apply the Bogoliubov method to our model. 
Within this approximation we first determine a conden- 
sate wave function and an excitation spectrum. Using 
the Bogoliubov spectrum we calculate occupation of the 
condensate and its fluctuation at finite temperatures. We 
compare these mean-field results to the results obtained 
within the exact model. We finish in Sec. IV with some 
concluding remarks. 
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II. EXACT RESULTS 

A. Ground state and excitation spectrum 

In our previous paper |lf| we have shown the algebraic 
method of diagonalization of the Hamiltonian describ- 
ing a system of many particles interacting via harmonic 
forces. The system under consideration consists of many 
particles confined by an external harmonic potential in- 
teracting by harmonic forces, i.e., two body interaction 
potential has the form: 

T/(x l -x,) = |r» 2 (x l -x J ) 2 , (1) 

where defines the interaction strength and a = +1 
signifies the attractive interaction of particles placed at 
positions x^ and x^- whereas a = — 1 - corresponds to 
repulsive interactions. The total Hamiltonian of the N- 
particle system has therefore the following form: 

N 1 

h =E2 (p ^ -+^)+t, v ^~^- ( 2 ) 

i=l i<j 

Let us first recall the exact results of . For the sake 
of simplicity we denote the set of all particle positions 
vectors by X^r = (xi, . . . , xjy). The Hamiltonian can be 
easily diagonalized if one introduces collective variables: 

X C N = Q^X N , (3) 

where X C N = (xj, . . . , x c N ) and the matrix Qjy = {<zjj } is 
orthogonal. One of these collective variables namely the 
center of mass of TV-particle system plays a particularly 
important role: 

1 N 
v i— l 

The choice of N — 1 remaining collective variables 
~X C N _ 1 — (x°, .. . ,x^ r _ 1 ) is not unique but this does not 
lead to any physical implications. In particular: 

N-l N 

(x^) 2 = E( x *) 2 = E x * - ( x ^) 2 - ( 5 ) 
i=i i=i 

In the following we are going to use a similar notation for 
description of a subsystem of s-particles, s = 1,...,N. 

The above-defined transformation brings the Hamil- 
tonian to the diagonal form and its eigcncncrgics can 
be easily found. While determining a spectrum, how- 
ever, one must take into account the proper symmetry 
of a total wave function. In the case of bosonic particles 
(N > 2) the allowed energies are: 

E=(~+m) + (?-(N-l)+n)u>, (6) 



where m = 0, 1, 2 . . ., n = 0, 2, 3 . . . and uj = y/l + aNti 2 . 
The first term describes excitations of the center of mass, 
i.e., d-dimensional harmonic oscillator of frequency equal 
to one. The second term in the Eq. (Q) corresponds to 
excitations of iV — 1 relative degrees of freedom. The fre- 
quency u> characterizes some effective potential felt by an 
individual quasi-particle because it results from a com- 
bined effect of all particles of our system. Let us observe 
that u> = 1 corresponds to the noninteracting case, the re- 
pulsive interactions give < lu < 1 while attractive forces 
lead to to > 1. Moreover, very small values of u> ss sig- 
nify very strong repulsion which almost destabilizes the 
whole system. It is very convenient to parameterize uj by 
an exponent k defined in the following way: 

lu = N K . (7) 

This exponent can be related to the actual strength of the 
interaction. In fact, for weakly interacting gas (a» « 1) we 
obtain very small values of this parameter: k ~ 0, while 
for strong interactions (cu « - repulsion, uj » 1 - at- 
traction) we have |«| 3> 1. Moreover, k is positive in the 
case of attraction while it is negative for repulsion. Let 
us add at this point that in realistic situations of short- 
range interparticle interactions, large Bose-Einstein con- 
densates can exist only for repulsive forces. In the case 
of attraction the size of the trapped condensate is lim- 
ited to about 1500 atoms g]. In our oscillatory model 
the forces between particles are negligible at small dis- 
tances, therefore the model leads to the condensation (in 
the thermodynamic limit) in both attractive and repul- 
sive case. 

The ground state of the system is the following: 

f (Xat) = *o(VwXw-i)*o(xSr), (8) 

where (X^ r _ 1 ,x^ r ) = <2jvX./v and the function 
$o(v / ^^jv-i) corresponds to the ground state of a sys- 
tem of N — 1 independent quasi-particles (in d spatial 
dimensions) interacting with an external potential of the 
harmonic oscillator of frequency w: 

$o(VCXSr_i) = (-) exp [-c(X$ V _ 1 ) 2 /2] , (9) 

and $o(xJ vr ) is the ground state of the single particle (cen- 
ter of mass) trapped into harmonic potential: 

*o(x$v)= (~) 1 exp [~(x^) 2 /2] . (10) 

Construction of excited eigenstates is difficult because it 
is not easy to impose the desired symmetry on the wave 
function. Such a procedure was describe in details in [^5| . 
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B. Order parameter and quantum depletion 

If the energy of the system (or equivalently the tem- 
perature) is sufficiently small we expect that the system 
forms a Bose-Einstein condensate. The BEC of the ideal 
gas manifests itself by a macroscopic occupation of the 
single particle ground state. In the case of interacting 
system it is not obvious what is this particular state 
which is 'macroscopically occupied'. The identification of 
the macroscopically occupied quantum state is equivalent 
to the definition of the order parameter - the single parti- 
cle wave function which is inherently related to the Bose 
condensation. The condensate subsystem can be then 
quite accurately described by the -/Vo-fold product of the 
order parameter, where No — O(N) is the occupation 
of the condensate. In the conventional approaches, for 
example in the Bogoliubov method, it is simply assumed 
that a mean value of the boson field operator is differ- 
ent than zero and this mean value is associated with the 
macroscopically occupied state. Then, consistently with 
the above assumption, the Bogoliubov equations give in 
fact the nonzero solution for the order parameter. How- 
ever, because of the superselection rules (resulting from 
the conservation of the barionic charge) any iV-particle 
system must be in the Fock state - the state with a well 
defined particle number. Therefore the mean value of the 
boson field operator must vanish in this state as the field 
operator changes the number of particles. 

In the following we use our model to demonstrate how 
to define the order parameter, occupation of the conden- 
sate, and its fluctuations. At zero temperature the sys- 
tem is in the ground state and one might naively expect 
that it is totally Bose condensated. However, the ground 
state of the N-particle bosonic system is not equivalent 
to the Bose-Einstein condensate. Interactions can signif- 
icantly deplete the condensate. We are going to show 
this effect in the most spectacular but also in relatively 
simple case of the zero temperature. 

Let us now define the hierarchy of the reduced s- 
particle density matrices which can be conventionally ob- 
tained by averaging the density matrix of the total system 
of N particles over the degrees of freedom of N — s re- 
maining particles. For a given iV-particle quantum state 
^(Xyv) the corresponding s-particle reduced density ma- 
trix p s (X s ; Y s ) is defined by: 



p s (X s ; Y s 



dR w _ s **(X s , Rjv_ s )*(Y s , R N s 



(11) 



We use previously defined shorthand notation for vectors 
in a configuration space of s-particles. The reduced den- 
sity matrix describes the subsystem of s-particles and can 
be directly related to different measurement processes. 
For the statistical description of the system one should 



first of all define the statistical density matrix by averag- 
ing all iV-particle density matrices with the appropriate 
statistical weights depending on the ensemble. In general 
it is quite a complicated task but at zero temperature 
there is only one quantum state of the system and no 
statistical averaging is necessary. 

The total wave function (or density matrix) carries all 
the information about the system. In real experiments 
however one does not probe simultaneously all the parti- 
cles. Typical detection scheme consists on the measure- 
ment of one or at most few particles at a given time. In 
other words a single measurement process is reduced to a 
subsystem of small number of particles. Such subsystems 
are described by reduced density matrices. In the con- 
sidered here case of zero temperature the description of 
the interacting system the s-particle density matrix can 
be brought to the following form: 



P*(X a ;Y s 



„CM 



(x^,y^)$ (v / ^X:_ 1 )cf ( v ^Y^ 1 ) 



(12) 



The functions <I>o describes the ground state of s — 1 
quasi-particles (collective relative coordinates, see Eq. 
(0)) while p CM corresponds to the density matrix of cen- 
ter of mass of the subsystem: 



„CM/ 



fUJ s \ d / 2 


S s c c 


w) exp 





exp 



-c\2 



c\21 



(yS) 



(13) 



The s-particles collective coordinates are defined in the 
familiar way: (X^.xJ) = Q S X S , (Y^y-) = Q S Y S 
and frequencies u> s , S s as well as auxiliary parameter "f s 
are: 



= 1 - 

UJ 

Is 



S(I-LU) 



N 



1 



N 



s(N 



7s 



(14) 
(15) 

(16) 



Having defined the s-particle matrices we are ready 
now to analyze the nature of the Bose-Einstein condensa- 
tion of the interacting system and to discuss the meaning 
of the order parameter. To this end we write the density 
matrix Eq. (n~2) in the diagonal form: 



Ps (X s ;Y s ) = ^AL s) ^ ) (X s )^ ) (Y s ). 



(17) 



The function 0^(X S ) can be treated as the wave func- 
tion of the s-particle subsystem: 



>(X,)=*o(VwXS-i)*»(>A£x£), 



(18) 



3 



where $o(V^-^s-i) ^ s tne ground state wave function 
of the relative degrees of freedom. This function cor- 
responds to the ground state of s — 1 noninteracting 
quasi-particles (in d-spatial dimensions) subject to the 
external harmonic potential of frequency u. The sec- 
ond part of the Eq. ( |l8| ) describes states of the center 
of mass of s-particles; $ n is simply the d-dimensional 
harmonic oscillator wave function corresponding to the 
effective center of mass frequency a s . Quantum numbers 
n = (ni, . . . ,nd) label different states of the center of 
mass while n = ni + . . . + rid corresponds to the energy 
of the given state. The effective center of mass frequency 



a s is: 



= [u> s (w s + 8 S )] 



1/2 



(19) 



It is interesting to observe that all the frequencies of the 
relative motion of the s-particles subsystem are exactly 
the same as the frequencies of the relative motion of the 
whole system, i.e., equal to u>. On the other hand the 
center of mass oscillation frequency of the subsystem is 
neither equal to lo nor to 1 (trap frequency). This col- 
lective degree of freedom couples to the center of mass 
of N — s remaining particles what leads to some energy 

(s) 

shift. Finally, the eigenvalues A n of p s are equal to the 
occupation probabilities of a given s-particle state: 



A« = 



d/2 /0 . 



(20) 



It follows from the normalization condition for the den- 
sity matrix that ^2 n An = 1- 

The spectral decomposition of the reduced single par- 
ticle density matrix gives natural single-particles states 
4>n (x). These states are crucial for the definition of the 
condensate wave function (order parameter). It can be 
seen from Eq.(poT) that if A goes to infinity (thermody- 
namic limit) with fixed value of the interaction frequency 
u) the lowest eigenvalue X^ dominates the others: 



A (1) 
A (1) 



AN 



if n + 0. 



(21) 
(22) 



This behavior signifies nothing else but the onset of the 
Bose-Einstein condensation. The single particle density 
matrix becomes very close to the pure state because 
with quite good accuracy it can be approximated by 
pi(x,y) « 



vXj 



q (y) . This particular single-particle 
■q (x) is usually called the order param- 



ground state 

eter. The TV-particle wave function can be quite accu- 
rately approximated by the iV-fold product of the order 
parameter. 

Above the critical temperature the situation is com- 
pletely different, namely all eigenvalues of pi should be 
close to zero what means that the single particle reduced 



density matrix is 'very far' from the pure state and or- 
der parameter vanishes - there is no leading state in the 
spectral decomposition of the single-particle density ma- 
trix. 




FIG. 1. Mean occupation of the condensate plotted as 
a function of the parameter k — log u/ log N for different 
number of particles; N = 10 3 - dashed line; N = 10 5 - dotted 
line, and N = 10 8 - full line. 

Our analytic formula allows to study quantitatively the 
role of interactions on the Bose-Einstein condensate. On 
the basis of the discussion it is obvious that the average 
occupation of the condensate becomes: 

(N ) = N [ dxdy4 1) (x)p 1 (x,y)4 1) (y) = NX^. (23) 



In the case of the ideal gas at zero temperature the above 
equation gives, of course, (No) = N; all particles occupy 
the single particle ground state. For a fixed number of 
particles, if the interaction strength | logo;] grows, the 
occupation of the condensate decreases. This behavior is 
presented in Fig. 1 where we show the mean occupation 
of the condensate versus the exponent k = log w/ log 
for different values of particle number A^ in three spa- 
tial dimensions (d — 3). The values of K less than zero 
signify repulsive interactions while K > corresponds 
to attraction. One can easily see that if the interaction 
becomes strong (\k\ ~ 1) the condensate is almost to- 
tally depleted. All curves presented in the figure tend to 
an universal curve if the number of particles increases. 
When N increases to infinity with k being constant then 
our expression for the occupation of the condensate has 
the form: 



(No) 
N 



i + ViV K - 1 + iv-( K + 1 ) + 1 



(24) 



The above formula, valid in the thermodynamic limit, 
gives an universal critical behavior. It exhibits no deple- 
tion (Ao = N) for | | < 1 followed by an abrupt jump 
and total destruction of the condensate (Nq = 0) for 
Id > 1. 
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The effect of quantum depletion of the trapped atomic 
condensate with a short range interactions, for the real- 
istic experimental parameters, has been estimated to be 
of the order of 1% j2|. This is opposite to the case of su- 
pcrfluid helium where this effect accounts for depletion 
as large as more than 90% [jl6| . Our model exhibits very 
interesting feature. It shows that in large N limit the 
quantum effects are almost negligible or totally destroy 
the condensate depending on the value of the interaction 
strength. At this point it is not clear if this is an unique 
feature of our model or if it is a more general result. 

We see that interactions play an important role. If 
they are strong, the condensate disappears although the 
TV-particle system remains in its ground state. There is 
no coherence in the strongly interacting system, i.e., no 
wave function can be assign to a single particle subsys- 
tem. At this point we want to make a comment about the 
notion of the coherence of the Bose-Einstein condensate. 
Approximate methods assume explicitly that the mean 
value of the boson field operator is different from zero in 
the case of the Bose-Einstein condensate. Therefore, the 
folk wisdom associates the condensate with the coherent 
state - the analog of the coherent state of the electromag- 
netic field. This analogy is of limited value and in fact 
may be misleading because the condensate must be in a 
Fock state in which a mean value the field operator van- 
ishes. However there is coherence in the condensate in the 
sense that majority of particles are described by the same 
wave function with the same phase. The expression for 
this wave function can be obtained rigorously only when 
one considers the single particle reduced density matrix. 

The 2-particle reduced density matrix allows to find 
a joint probability of finding one particle in a given sin- 
gle particle state and simultaneously another particle in 
another given state. In particular we have: 

(N (N -1)) = N(N-1) y"dX 2 dY 2j02 (X 2 ,Y 2 ) 

^ 1) (x 1 )^ 1) (x 2 )^ 1) (y 2 )4 1) (y 1 ). (25) 
Simple integration gives: 

M-D)=iv(iv-i)(^) d (^) d 

\LJ + Oil J \&2 + Oi\ J 
/ i \ d / 2 

x — T ~r~F~ • ( 26 ) 

\w 2 + cti + d 2 J 

Now we are ready to analyze the fluctuations of the con- 
densate defined as: 

(6 2 N ) = (A 2 ) - <A ) 2 . (27) 

These fluctuations are shown in Fig. 2. We see that as 
the interaction strength grows up (at fixed number of 
particles) the fluctuations start to grow from zero value 
for the ideal gas. However when the interactions become 



so strong that condensate practically disappears (\k\ ~ 1) 
fluctuations also decrease - as there is no condensate the 
fluctuations also die out. The fluctuations are maximal 
in a region of the critical destruction of the condensate 
by quantum effects. 
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FIG. 2. Fluctuations of the condensate plotted as a func- 
tion of the parameter « = log a;/ log iV for different number 
of particles; N = 10 3 - dashed line; N = 10 5 - dotted line, 
and N = 10 8 - full line. 

C. Finite temperatures 

In this subsection we estimate, within the exact model, 
some effects in finite-temperature behavior of the trapped 
gas. Rigorous description of the condensate requires a 
knowledge of the statistical density matrix of the N- 
particle system. Knowing this matrix one can apply 
the procedure described previously to define the finite 
temperature condensate, its occupation and fluctuations. 
However, because of a huge degeneracy of high energy 
states the statistical averaging procedure is difficult. 
Therefore, we limit our study to the case of weak interac- 
tions < 1) when we can neglect the quantum effects 
leading to a significant depletion of the condensate. In 
this case we can expect that the condensate wave func- 
tion in a finite temperature is equal to the ground state 
of the harmonic oscillator with some effective frequency 
(characterizing a mean field experienced by a single par- 
ticle) which at T — is equal to a%. In the case of a 
weak interaction this frequency can be approximated by 
ax ~ 10 ■ We expect that, similarly to the noninteract- 
ing case, the main effect of the temperature is to deplete 
the condensate rather than modify the condensate wave 
function, i.e., the frequency ol\. 

As it has been shown in jL5| , the trace of the density 
matrix of the A-particle system, i.e., the microcanoni- 
cal partition function T(N,E) is identical (in the ther- 
modynamic limit) to the microcanonical partition func- 
tion Tq(N, E,cu) of the N noninteracting bosonic parti- 
cles (quasi particles) trapped by the harmonic potential 
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of frequency u: 



T(N,E) w T (N,E,w). 



(28) 



The effect of the center-of-mass excitations on the spec- 
trum and on the statistical properties of the system is 
negligible since it is related to only one degree of free- 
dom as compared to the N — 1 remaining collective de- 
grees of freedom. This fact, together with our remarks 
about the condensate wave function, signifies that the 
system of interacting (via harmonic forces) particles is, 
in the thermodynamic limit, equivalent to the ideal gas 
in the oscillatory trap. 

This observation allows to recall all results obtained 
for the ideal Bose gas p7| . In particular, considered here 
interacting system undergoes the Bose-Einstein conden- 
sation at the temperature T c equal to: 



N 

W) 



1/3 



(29) 



where the ( is the Riemann function, £(3) = 1.2020569. 
This critical temperature should be compared to the crit- 
ical temperature of the noninteracting system, 



Tn = 



( N 



VC(3) 



1/3 



The shift of the critical temperature AT = T c 
therefore equal to: 



AT 



= oj-1. 



(30) 



T is 



(31) 



It is negative in the case of repulsive interaction and has 
the opposite sign for the attractive system what is in a 
qualitative agreement with results of Ref. (ll| for the sys- 
tem of trapped atoms. The critical temperature for the 
trapped gas with repulsive interactions is decreased as 
compared to the noninteracting case: due to interactions 
a mean separation between particles grows, therefore the 
quantum statistical effects become important at larger 
de Broglie wavelength. 

It is worth to stress that the transition temperature 
for the system with short range interactions remains a 
controversial subject even in the uniform case p9| . There 
is no consensus how the shift of the temperature should 
depend on the interaction strength, nor even the sign. 

The fraction of condensate particles (No) is: 



The value of the fluctuations of the number of the par- 
ticles in the condensate calculated for different statis- 
tical ensembles differ significantly ^0|. In the canoni- 
cal ensemble the square of the fluctuations (<5 2 A*o)cn = 



CN 



{N 



,2 
CN 



is: 



(<5 2 iV )cN = 



T 



C(2), 



(33) 



where £(2) = 7r 2 /6. In the case of perfectly isolated 
system (microcanonical ensemble) fluctuations of the 
gro und state occupation are smaller than canonical ones 

C«b>» -(§)'«« (l-j^j). m 
and C(4) = tt 4 /90. 



III. BOGOLIUBOV APPROXIMATION 

A. Ground state 

In realistic cases the exact solution of the many-body 
Schrodinger equation is impossible. Instead, mean-field 
approaches are being developed. The basic idea of the 
mean-field theory was introduced by Bogoliubov. We will 
now formulate the Bogoliubov approximation in the case 
of inhomogeneous systems [2lJ . Next we find the energy 
spectrum within this approximation for our exactly solv- 
able model. In order to do this we rewrite the Hamilto- 
nian (||) using the second quantization formalism. Thus 
we introduce the field operator ip(x) which annihilates 
a particle at a point x and its conjugate V>^(x) which 
creates a particle at a point x. These operators fulfill 
standard bosonic commutation relations: 



ip(x),ft(x') =(J(x-x'). 



(35) 



In the second quantization formalism the Hamiltonian 
becomes: 



H = J dx^ t (x)iJ (x)?A(x) 



-- / dxdx> t (x)v3 t (x')^(x-x')^(x)Vi(x'), (36) 



where 



(32) 



It is known that the critical temperature and mean oc- 
cupation of the ground state do not depend on the sta- 
tistical ensemble used for the description of the system. 



flo(x) 



;(p 2 + x 2 ). 



(37) 



The Bogoliubov approximation is formulated in two 
steps. The first one is to express the field operator as a 
sum of its mean value y / A'o0o( x ) & n d an operator </>(x) 
responsible for the fluctuations around the mean value. 
Below the Bose-Einstein condensation temperature the 
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occupation of the ground state is nonzero. It is therefore 
convenient to write 



4 1} (v^x) 



ai\z ,1 2 
exp(--aix ), 



(43) 



^(x) = v/iVo0o(x) + (/>(x). 



(38) 



where the effective frequency a\ in the limit of weak in- 
teractions can be approximated by: 



Note that now the mean value of the field operator ^(x) 
is not equal to zero and is proportional to the square root 
of the number of condensed particles No- The spirit of 
the Bogoliubov approximation is based on the fact that 
the occupation of the condensate is of the order of to- 
tal particles number No ~ O(N). This, in principle, 
limits the Bogoliubov approach to low temperatures. In 
the following, consistently with the above assumption we 
will substitute in the Eq. (|38|) Nq by N. For any phys- 
ically realizable ./V-particlc state, and in particular for 
our solutions of the system Eq. (||) , the mean value of the 
field operator is zero. (Strictly, separation of the Eq.(|3"§|) 
should be done for operators conserving particle number. 
However, extracting a c-number part of the field oper- 
ator has formally the same consequences and is easier 
to handle). If such a form is substituted into the total 
Hamiltonian Eq. ([36]) the number of particles is no longer 
conserved. To overcome this difficulty one considers the 
grand canonical Hamiltonian instead 



cii 1 



1 - U) 
N 



(44) 



K = H — fj,N, 



(39) 



where N is the total particle number operator and fi is 
the chemical potential. It should be chosen in such a way 
that the mean particle number is equal to the desired 
value. Self consistent equation for the condensate wave 
function 0o( x ) follows from the assumption that the de- 
composition ( |38| ) gives the best self consistent function. 
We find: 

{ff o (x) + K//[0o,x]}0o(x) =M)(x), (40) 
where the effective potential is: 



Because the effective frequency a± is very close to lj the 
Bogoliubov expression for the condensate wave function 
is quite accurate in the limit of weak interactions. How- 
ever, if the interaction strength is large the Bogoliubov 
method fails to reproduce the condensate wave func- 
tion. This is consistent with the basic assumption of the 
Bogoliubov method which requires the order parameter 
(multiplied by the mean occupation of the condensate) 
to be large. As we have shown in the previous section 
this is not the case for strongly interacting system. 



B. Excitation spectra 

The second step in the Bogoliubov method is to find 
the low energy excitation spectrum of the system by ex- 
panding the total Hamiltonian around the mean value 
of the field operator given by the solution Eq.(42). Af- 
ter substituting the field operator Eq.(|3^) into K and 
retaining all terms up to 0((j) 2 ), the operator K can be 
diagonalized with the help of a canonical transformation: 



A 



u x (x)(3 x + v* x (x)p{ 



(45) 



where (3\ and j3\ are bosonic annihilation and creation 
operators. The diagonal form of the operator K is 

K = ^A X f dxuX(x)« A (x) + Y, A a/?!/3a, (46) 

A A 

V ef f[(j>o,x] = T^NCi 2 J dxVS(x')(x - x') 2 ^o(x')- (41) provided that functions C/ A (x) and V x (x) defined as: 



This equation replaces the standard Gross-Pitaevskii 
equation for the condensate wave function. The effective 
potential in the Eq.(^o|) has different form than usuall 
nonlinear term appearing in the Gross-Pitaevskii equa- 
tion because of long range forces assumed in our model 
as opposed to the more realistic zero range interactions. 

The lowest energy state of the Hamiltonian from 
Eq.© is 



U\(x) = u\(x) + u A (x), 
V\(x) = it A (x) - u A (x), 

satisfy the normal-mode equations: 



(47) 
(48) 



<Mx) = (-) 4 exp(-icjx 2 ), 



(42) 



and the value of chemical potential is fi = 3/4(l+a>). The 
function Eq. ( 42 ) is the Bogoliubov approximation to the 
exact order parameter <j)\p (-^/aJx) found in the previous 
section: 



H u {x)U x (x) + J dx'G(x,x')C/ A (x') - A A y A (x), (49) 
H w (x)Vx(x) = A A C/ A (x), (50) 
with the following normalization condition: 

dx(?4(x)M A /(x) -vl(x)vx'(x)) = 6x,\>. (51) 

In these formulas A A has a meaning of an eigenvalue, 
H w (x) — l/2(p 2 +u> 2 x 2 ) — 3lo/2, and the integral kernel 
G(x,x') is: 
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G(x, x') = 2JV<^,(x)y(x - x')<M*')- 



(52) 



The above equations can be easily solved if we expand 
functions U\(x.) and Vx(x) in the basis of eigenfunctions 
V> n (x) = ipnxnynAx-) of the Hamiltonian i? w (x): 



[7 A (x) = ^a^n(V^x), 



F A (X)= Vfe^n(V^x), 



(53) 
(54) 



where components of the vector n = (n x ,n y ,n z ) are 
the standard quantum numbers of the harmonic oscil- 
lator eigenfunction of energy equal to ruo, where n = 



The coefficients and have to be determined from 
the Eqs. ( p9|) and (|50|). Let us remind that the solution 
of the Eq. (|40|) for the order parameter </>o(x) is the first 
function of the chosen basis set, <^o(x) = V'ooo(\/^' x )- 
For this reason and also due to the oscillatory form of 
the inter-particle interactions, the integral kernel in Eqs. 
( [49| ) and ( |50| ) couples only these basis functions tp n which 
correspond to three lowest eigenstates (n = 0, 1 and 2) 
of the Hamiltonian iJ w (x). For larger n quasi-particles 
excitation energies are those of the harmonic spectrum: 



A„ = 



nuj, 



if n > 2. 



(55) 



These eigenvalues are degenerated. In general, there is a 
close link between the number of eigenmodes «a( x ) and 
v\ (x) corresponding to the eigenvalue nu> and the number 
of the oscillatory states of the same energy. Therefore, in 
order to classify independent solutions, it is convenient to 
use oscillatory quantum numbers instead the parameter 
A which simply enumerates quasi-particles eigenmodes. 
With this notational modification the solutions of the 
Eq. ( (49| , [50|) corresponding to energies A„ with n > 2 are 
the following: 



U n (x) = Vnlv^x), 
Vn(x) = 0. 



(56) 
(57) 



The low laying interacting states involve coupling of the 
bare oscillatory eigenfunctions. For the lowest excitation 
energy we get: 



A 1 = l. 



(58) 



There are three different modes of that energy corre- 
sponding to the excitation of one of the x, ?/, or z de- 
gree of freedom. Below we present only one pair of the 
eigenmodes as the remaining two can be obtained by the 
permutation of the indices only. The x-direction eigen- 
modes are: 



Mioo(x) = — 7 =^Vioo(V^x), 



"100 (x) 



'4cj 



■0100 (Vw> 



(59) 
(60) 



Let us notice that energy of this mode of excitations is 
equal to the single excitation quantum of the trap mode. 

The second excitation energy is equal to one of the 
excitation energy of the Hamiltonian 7J w (x), namely: 



2u. 



(61) 



There are six different eigenmodes of the above energy 
which is exactly the degeneracy of the second state of 
the 3D oscillator. The first three pairs of them are re- 
lated to the double excitation along one of the axis of the 
coordinate system and are of the form: 

_ UJ 2 -1 _ 

M2oo(x) = , 02oo(v / ^x) + j= Vooo(v^x), (62) 
2v 2ui 



«200(x) 



oj z - 1 



2V2< 



^ooo(V^x), 



(63) 



(the other two pairs of eigenstates can obtained by the 
permutation of the oscillatory quantum numbers, as pre- 
viously). The remaining three pairs of eigenmodes cor- 
respond to two single quanta of excitations along two 
different principal axis of the coordinate system. For ex- 
ample, one such pair is: 



itno(x) = ipuo(^/ux), 
viio(x) = 0, 



(64) 
(65) 



and two others can be obtained by permutations of in- 
dices. 

The shift in the ground state energy is given by 
(1-c) 2 



^ A n J dxv 2 x (x) = 



4w 



-(2 + 2w + uj 2 ). (66) 



The excitation energies A„ obtained within the Bo- 
goliubov approach are the same as exact eigenenergies of 
the interacting system. The degeneracies of the eigenen- 
ergy state when we compare with Ref. Jl5| are also the 
same. The Bogoliubov method is very well suited for 
the description of the excitation spectrum of the quan- 
tum degenerate gas. This result is somewhat surpris- 
ing. One might rather expect that Bogoliubov approach 
works well only for the short range interaction. Our cal- 
culation shows that it works also in a rather exotic case 
when the interaction strength grows quadratically with 
the distance between particles. There are some differ- 
ences between wave functions obtained in the Bogoli- 
ubov approximation and exact solutions of the TV-particle 
Hamiltonian They come from the fact that the ex- 
act solution cannot be written as a symmetrized product 
of any single particle functions. 



C. Condensate fraction and fluctuations 

In this subsection we consider Bose gas at finite tem- 
peratures. We will study the impact of interactions on 
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the occupation of the condensate and its fluctuations us- 
ing Bogoliubov method p2fl . We will compare the ob- 
tained results with those inferred from the exactly soluble 
model. 

The statistical density matrix is p = Z^ 1 exp(—K/T), 
where K is the Hamiltonian given in Eq.(|3|) and Z is 
the statistical sum. This density matrix describes the ex- 
cited subsystem only (quasi-particles) . Number of quasi- 
particles is not conserved and the condensate is assumed 
to act as a reservoir of quasi-particles. In fact all this 
assumptions are in the spirit of the Maxwell's demon en- 
semble introduced for the description of the ideal gas 
below the condensation temperature (20|. Imposing the 
constraint on the total number of particles implies that 
occupation and fluctuations of the condensate can be di- 
rectly related to the mean number and fluctuations of 
quasi-particles. 

The mean number of excitations (quasi particles) 
above the condensed phase (N e ) = N — (No) is defined 
as: 

(N e ) = J dx(£t(x)£(x)), (67) 
what leads to the following expression: 



{N e 



E 



J dx { [ul (x) + v 2 n (x)] in + v 2 n (x) } , (68) 



where /„ = [exp(A„/T) — plays the role of the mean 
quasi-particle occupation of the given energy state. The 
functions u n and v n are closely related to the wave func- 
tions of harmonic oscillator. The integration can be easily 
performed but the summation over all eigenstates might 
be quite difficult because of huge degeneracy of the en- 
ergy levels. Fortunately, in our case almost all functions 
i> n vanish and their contribution to the final result is neg- 
ligible at finite temperatures. Therefore the problem can 
be easily reduced to the calculation of the canonical oc- 
cupation of the condensate trapped in the harmonic trap 
of frequency u>. Again, similarly as in the exact solu- 
tion we will treat separately the two regimes: (i) zero 
temperature limit where only quantum effects described 
by the last term of the Eq. ( p8[ ) affect the condensate 
population, (ii) finite temperature case where the above 
mentioned effect is negligible. The occupation of the in- 
teracting condensate at zero temperature is: 

(No) =N-^( 1 - ^ - 8^a(l - ^) 2 - (69) 

The first term in Eq. (|69| ) corresponds to the weak inter- 
action limit of the exact result, which is 



(No) 



(70) 



The second term in Eq. (^9|) is new. This additional 
term overestimates the depletion of the condensate as 



compared to the rigorous treatment. Therefore the Bo- 
goliubov treatment does not describe quantitatively an 
occupation of the zero temperature condensate. 

On the other hand, in the thermodynamic limit, we 
recover the familiar expression: 



(N 



N 



T 



C(3). 



(71) 



This is exactly the same result which we obtained in 
the exact treatment. Similarly, the critical temperature 
which can be defined by setting (No) to zero in the Eq. 
7 H): 



T r 



t N \ 1/3 
U(3)J 



(72) 



is identical to the exact result, Eq. (|29|). It might be 
somewhat surprising but the Bogoliubov method works 
in our case pretty well up to the critical temperature, i.e., 
in the region where, in principle, its assumptions are not 
valid. 

In a similar way we can obtain the fluctuations of the 
condensate population. Because fluctuations, contrary to 
the mean occupation, depend on the statistical ensemble, 
we use a notation which explicitly indicates the kind of 
performed averages (i.e., the canonical one). As we im- 
posed constraint on the total number of particles, these 
fluctuations are equal to the fluctuations of the above- 



condensate part, (SN§)cn = (Nq 



CN 



CN 



W/dx [<x n (x)+« n (x)] 2 (/„• + /„) 



dxu n (x)u n (x) 



(/n + /n + 1) >, (73) 



where the last term is responsible for the quantum fluc- 
tuations which do not vanish at zero temperature. Again 
the problem of calculating the condensate fluctuation can 
be reduced to finding the fluctuations of the ideal trapped 
Bose gas. In the case of zero temperature we have: 



(<5iV 2 )cN 



2^ 



2 (1 



,2\2 



1 



(74) 



This formula overestimates (quite significantly) the con- 
densate fluctuations due to the quantum effects as com- 
pared to the exact results in the limit of weak interaction, 
Eq. @. 

The fluctuations (in the thermodynamic limit and af- 
ter neglecting the small quantum fluctuations obtained 
above) at finite temperatures are: 



tf.Y,?)rN = (~) ^<->- 



(75) 



This result gives the correct value of the thermal fluctu- 
ations of the condensate. The Bogoliubov method works 
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very well in predicting the thermal fluctuations of the 
condensate while it fails to reproduce the correct value 
of the quantum fluctuation. It might be surprising as in 
fact the approach should work in the low temperature re- 
gion. In our opinion, it is the substitution of the operator 
annihilating the lowest energy state by a c-number that 
is responsible for inaccurate treatment of some quantum 
effects, particularly important at zero temperature. 

IV. CONCLUSIONS 

In our paper we used the exactly soluble many-particle 
model to illustrate the rigorous procedure of defining the 
condensate phase at zero temperature. By diagonaliz- 
ing one-particle reduced density matrix we were able to 
study in details the role of interactions on the condensate 
at zero temperature. If the interaction strength becomes 
large |k| > 1 condensate disappears even when the sys- 
tem is in its ground state. This total depletion of the 
condensate has, in the thermodynamic limit, a character 
of critical phenomena. The destruction of the condensate 
signifies a breakdown of the mean-field theory. Due to 
strong quantum correlations the system cannot be viewed 
as being composed of independent quasi-particles mov- 
ing in some effective potential resulting from interactions 
with a rest of the system. Instead the quasi-particles be- 
come strongly correlated and simple single-particle pic- 
ture is not longer valid. 

We have also carefully compared the exact quantum 
solutions of the oscillatory model with the approximate 
solutions obtained with the help of the Gross-Pitaevskii 
equation and Bogoliubov approximation. We have found 
that many of the characteristics of the exact solutions 
like the excitation spectrum, occupation of the conden- 
sate, and its thermal fluctuations are indeed reproduced 
with the help of the approximate methods. The Bogoli- 
ubov approach fails in the case of very strong interactions 
when the condensate is almost destroyed. This is how- 
ever consistent with the basic assumption of the Bogoli- 
ubov method which explicitly assumes small condensate 
depletion and relays on the validity of a mean-field de- 
scription. Surprisingly, in the studied model, the method 
works quite well even at temperatures close to the criti- 
cal one. On the other hand the zero temperature (quan- 
tum) depletion and fluctuations of the condensate are not 
given correctly by the Bogoliubov method. Substitution 
of the destruction operator of the particle in the low- 
est energy state by a classical field results in inaccurate 
description of some quantum effects. Fortunately for the 
weakly interacting system, these effects are small and can 
be neglected. 
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